Gam

## Loaded package itsadug 1.0.2 (see 'help("itsadug")' ).

## 
## Attaching package: 'itsadug'

## The following object is masked from 'package:ggplot2':
## 
##     alpha

e <- read_excel("~/SharedFiles/ST606/2020/data/Exercise/fit_database_anthropometric_all.xlsx")

#Dealing with NAs
ena<-na.omit(e)

#Renaming Variables(Renaming Column Names for easy code handling and Cleaning Data)
names(ena)[names(ena)=='measurement date']<- 'measurement_date'
names(ena)[names(ena)=='age (years)']<- 'age_years'
names(ena)[names(ena)=='age bin']<- 'age_bin'
names(ena)[length(ena)]<-"z_cat_WHO"
names(ena)[9]<-"z_score_WHO"
names(ena)[6]<-"height_cm"
names(ena)[7]<-"weight_kg"

#Removing "NA" Characters
ena$observation <- 1:nrow(ena)
x<-ena[ena$z_cat_WHO=="NA",]
y<-x$observation
ena<-ena[-y,]

#Adding additional column
ena$year <- year(ena$measurement_date)
ena$ID <- as.factor(ena$ID)
#Changing Data Types
ena <- ena %>%
  mutate(gender = as.factor(gender),
         z_cat_WHO=as.factor(z_cat_WHO),
         measurement_date=as.Date(measurement_date),
         BMI = as.numeric(BMI),
         z_score_WHO=as.numeric(z_score_WHO),
         height_cm=as.numeric(height_cm),
         weight_kg=as.numeric(weight_kg))

Diving data as per Year

ena_2007 <- filter(ena, year(measurement_date) == 2007)
ena_2008 <- filter(ena, year(measurement_date) == 2008)
ena_2009 <- filter(ena, year(measurement_date) == 2009)
ena_2010 <- filter(ena, year(measurement_date) == 2010)
ena_2011 <- filter(ena, year(measurement_date) == 2011)
ena_2012 <- filter(ena, year(measurement_date) == 2012)
ena_2013 <- filter(ena, year(measurement_date) == 2013)
ena_2014 <- filter(ena, year(measurement_date) == 2014)
ena_2015 <- filter(ena, year(measurement_date) == 2015)
ena_2016 <- filter(ena, year(measurement_date) == 2016)
ena_2017 <- filter(ena, year(measurement_date) == 2017)
ena_2018 <- filter(ena, year(measurement_date) == 2018)

Checking uniques IDs whose observation has been taken continuously from 2007 to 2018

#kids data of 2007 whose obeservation were taken till the end of 2018 without any miss.
a_07_08 <- subset(ena_2007, ID %in% ena_2008$ID)
a_07_09 <- subset(a_07_08, ID %in% ena_2009$ID)
a_07_10 <- subset(a_07_09, ID %in% ena_2010$ID)
a_07_11 <- subset(a_07_10, ID %in% ena_2011$ID)
a_07_12 <- subset(a_07_11, ID %in% ena_2012$ID)
a_07_13 <- subset(a_07_12, ID %in% ena_2013$ID)
a_07_14 <- subset(a_07_13, ID %in% ena_2014$ID)
a_07_15 <- subset(a_07_14, ID %in% ena_2015$ID)
a_07_16 <- subset(a_07_15, ID %in% ena_2016$ID)
a_07_17 <- subset(a_07_16, ID %in% ena_2017$ID)
a_07_18 <- subset(a_07_17, ID %in% ena_2018$ID)

#All data of the a_07_18 dataset kids from 2007 to 2018.
a_unique_07<- subset(ena, ID %in% a_07_18$ID)

#Checking prop increase per year
a_unique_07_pct<-a_unique_07 %>%
  group_by(ID) %>% 
  dplyr::arrange(measurement_date, .by_group = TRUE) %>%
  dplyr::mutate(pct_change_ht = (height_cm/lag(height_cm) - 1) * 100) %>% 
  dplyr::mutate(pct_change_wt = (weight_kg/lag(weight_kg) - 1) * 100) %>%
  dplyr::mutate(pct_change_bmi = (BMI/lag(BMI) - 1) * 100) %>%
  dplyr::select(ID, measurement_date, height_cm, BMI, weight_kg,pct_change_ht,pct_change_wt,pct_change_bmi,gender)

#Height
shortest_07_boy <- a_07_18 %>% arrange(height_cm) %>% filter((gender) %in% c("boy"))%>% head(5)
shortest_07_boy_IDs<-shortest_07_boy$ID

shortest_07_girl <- a_07_18 %>% arrange(height_cm) %>% filter((gender) %in% c("girl")) %>% head(5)
shortest_07_girl_IDs<-shortest_07_girl$ID

#Weight
weightl_07_boy <- a_07_18 %>% arrange(weight_kg) %>% filter((gender) %in% c("boy"))%>% head(5)
weightl_07_boy_IDs<-weightl_07_boy$ID

weightl_07_girl <- a_07_18 %>% arrange(weight_kg) %>% filter((gender) %in% c("girl")) %>% head(5)
weightl_07_girl_IDs<-weightl_07_girl$ID

#BMI
bmil_07_boy <- a_07_18 %>% arrange(BMI) %>% filter((gender) %in% c("boy"))%>% head(5)
bmil_07_boy_IDs<-bmil_07_boy$ID

bmil_07_girl <- a_07_18 %>% arrange(BMI) %>% filter((gender) %in% c("girl")) %>% head(5)
bmil_07_girl_IDs<-bmil_07_girl$ID

#Checking weight rate of the lightest kids from 2007 throught 2018
#BOY
lightest_boy<-a_unique_07_pct %>% filter((ID) %in% c(weightl_07_boy_IDs)) %>% 
  select(ID, measurement_date, weight_kg, pct_change_wt)

ggplot(lightest_boy, aes(x=measurement_date, y = pct_change_wt ))+
  geom_line()+facet_wrap(~ID)+
  geom_hline(yintercept=0, linetype="dashed", color = "red")

## Warning: Removed 1 rows containing missing values (geom_path).

#GIRL
lightest_girl<-a_unique_07_pct %>% filter((ID) %in% c(weightl_07_girl_IDs)) %>% 
  select(ID, measurement_date, weight_kg, pct_change_wt)

ggplot(lightest_girl, aes(x=measurement_date, y = pct_change_wt ))+
  geom_line()+facet_wrap(~ID)+
  geom_hline(yintercept=0, linetype="dashed", color = "red")

## Warning: Removed 1 rows containing missing values (geom_path).

#Checking growth rate of the shortest kids from 2007 throught 2018
#BOY


shortest_boy<-a_unique_07_pct %>% filter((ID) %in% c(shortest_07_boy_IDs)) %>% 
  select(ID, measurement_date, height_cm, pct_change_ht)

ggplot(shortest_boy, aes(x=measurement_date, y = pct_change_ht ))+
  geom_line()+facet_wrap(~ID)+
  geom_hline(yintercept=0, linetype="dashed", color = "red")

## Warning: Removed 1 rows containing missing values (geom_path).

#GIRL
shortest_girl<-a_unique_07_pct %>% filter((ID) %in% c(shortest_07_girl_IDs)) %>% 
  select(ID, measurement_date, height_cm, pct_change_ht)

ggplot(shortest_girl, aes(x=measurement_date, y = pct_change_ht ))+
  geom_line()+facet_wrap(~ID)+
  geom_hline(yintercept=0, linetype="dashed", color = "red")

## Warning: Removed 1 rows containing missing values (geom_path).

#Negative growth rate check (Wrong Data)
shortest_girl%>% filter((ID)%in%c(262))

## # A tibble: 21 x 4
## # Groups:   ID [13,540]
##    ID    measurement_date height_cm pct_change_ht
##    <fct> <date>               <dbl>         <dbl>
##  1 262   2007-10-01             116        NA    
##  2 262   2008-04-01             120         3.45 
##  3 262   2008-10-01             122         1.67 
##  4 262   2009-04-01             125         2.46 
##  5 262   2009-10-01             127         1.6  
##  6 262   2010-04-01             131         3.15 
##  7 262   2010-10-01             135         3.05 
##  8 262   2011-04-01             137         1.48 
##  9 262   2011-10-01             138         0.730
## 10 262   2012-04-01             141         2.17 
## # … with 11 more rows

Weight Analysis

a_unique_07$is_boy <- as.factor(ifelse(a_unique_07$gender=='boy', "1", "0"))

#Boys Weight
a_unique_07_boys <- filter(a_unique_07, gender=="boy")
ggplot(a_unique_07_boys, aes(x=age_years, y=weight_kg, color=ID)) + geom_line()

#Girls Weight
a_unique_07_girls <- filter(a_unique_07, gender=="girl")
ggplot(a_unique_07_girls, aes(x=age_years, y=weight_kg, color=ID)) + geom_line()

gam() or bam()

There are two functions for implementing a GAMM model: gam() and bam(). There are largely similar. The most important difference is that bam() is optimized for big data sets. We chose the Gam().

Smooths terms To model a potentially nonlinear smooth or surface, three different smooth functions are available: * s() : for modeling a 1-dimensional smooth, or for modeling isotropic interactions (variables are measured in same units and on same scale)

te(): for modeling 2- or n-dimensional interaction surfaces of variables that are not isotropic (but see info about d parameter below). Includes ‘main’ effects.

ti(): for modeling 2- or n-dimensional interaction surfaces that do not include the ‘main effects’.

Parameters of smooth functions The smooth functions have several parameters that could be set to change their behavior. The most often-used parameters are listed here:

k : number of ‘knots’. This parameter determines the upper bound of the number of underlying base functions being used to build up the curve. Thus, this parameter constraints the wigglyness of a smooth, or - as a metaphor - the number of bowpoints of a curve. Note that the model will base the number of base functions (reflected in the edf of the summary) on the data with the setting for k as upper bound. By default, the value of k for s() is around 9, and for te() and ti() 5 per dimension. Importantly, the value of k should be at most one less than the number of unique data points, otherwise it will fit the density of that predictor.

d : for specifiying that predictors in the interaction are on the same scale or dimension (only used in te() and ti()). For example, in te(Time, width, height, d=c(1,2)), with width and height reflecting the picture size measured in pixels, we specify that Time is on a different dimension than the next two variables. By default, the value would be d=c(1,1,1) in this case.

bs: specifies the type of underlying base functios. For s() this defaults to “tp” (thin plate regression spline) and for te() and ti() this defaults to “cr” (cubic regression spline). For random intercepts and linear random slopes use bs=“re”, but for random smooths use bs=“fs” (see below).

Setting up a GAMM model Before the different components of a GAMM model were summarized. In this section we focus on finding the model that best fits the data.

We start with the full model. Note that we did not include interactions with the predictor Trial, because it is not the focus of our simulated experiment, and it takes much longer to run. Therefore, we only included Trial as random effect. We did not include Item effects, because these were not specified in this data set.

Random effects Three different types of random effects are distinghuished when using GAMMs:

random intercepts adjust the height of other modelterms with a constant value: s(ID, bs=“re”).

random slopes adjust the slope ofthe trend of a numeric predictor: s(ID, age, bs=“re”).

random smooths adjust the trend of a numeric predictor in a nonlinear way: s(age, ID, bs=“fs”, m=1).

Notes:

Random intercepts and random slopes could be combined, but the random smooths already include random intercepts and random slope effects.

The argument m=1 sets a heavier penalty for the smooth moving away from 0, causing shrinkage to the mean.

No Random Effect

model1 <- bam(weight_kg ~ is_boy + age_years,
          data=a_unique_07)
summary(model1)

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## weight_kg ~ is_boy + age_years
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -5.04225    0.91078  -5.536 3.36e-08 ***
## is_boy1      2.61234    0.43921   5.948 3.03e-09 ***
## age_years    4.22479    0.06788  62.239  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## R-sq.(adj) =  0.567   Deviance explained = 56.7%
## -REML =  11718  Scale est. = 143.8     n = 3002

gam.check(model1)

## 
## Method: REML   Optimizer: outer newton
## full convergence after 5 iterations.
## Gradient range [-0.0010357,-0.0010357]
## (score 11718.01 & scale 143.7979).
## Hessian positive definite, eigenvalue range [1499.501,1499.501].
## Model rank =  3 / 3

Random intercepts Effect

model2 <- bam(weight_kg ~ is_boy + age_years
          + s(ID, bs="re"),
          data=a_unique_07)
#summary(m2)

plot(model2)

gam.check(model2)

## 
## Method: fREML   Optimizer: perf newton
## full convergence after 6 iterations.
## Gradient range [-3.708556e-05,2.774985e-05]
## (score 9940.748 & scale 36.21685).
## Hessian positive definite, eigenvalue range [64.11455,1502.763].
## Model rank =  144 / 144 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##        k' edf k-index p-value
## s(ID) 141 137      NA      NA

Random intercepts + slopes Effect

model3 <- bam(weight_kg ~ is_boy + age_years
          + s(ID, bs="re")
          + s(ID, age_years, bs="re"),
          data=a_unique_07)
#summary(m2)

plot(model3)

gam.check(model3)

## 
## Method: fREML   Optimizer: perf newton
## full convergence after 6 iterations.
## Gradient range [-3.222143e-06,2.737994e-06]
## (score 9175.449 & scale 18.11193).
## Hessian positive definite, eigenvalue range [48.11577,1505.282].
## Model rank =  285 / 285 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##                  k' edf k-index p-value
## s(ID)           141 125      NA      NA
## s(ID,age_years) 141 134      NA      NA

Comparing Models

library(itsadug)
compareML(model2,model3)

## model2: weight_kg ~ is_boy + age_years + s(ID, bs = "re")
## 
## model3: weight_kg ~ is_boy + age_years + s(ID, bs = "re") + s(ID, age_years, 
##  model2: weight_kg ~ is_boy + age_years + s(ID, bs = "re")
## 
## model3:     bs = "re")
## 
## Chi-square test of fREML scores
## -----
##        Model    Score Edf   Chisq    Df  p.value Sig.
##       model2 9940.748   4                            
## fREML model3 9175.449   5 765.299 1.000  < 2e-16  ***
## 
## AIC difference: 1966.43, model model3 has lower AIC.

Random intercepts + slopes + smooths Effect

model4 <- bam(weight_kg ~ is_boy + age_years
          + s(ID, bs="re")
          + s(ID, age_years, bs="re")
          + s(age_years, ID, bs="fs", m=1),
          data=a_unique_07)
#summary(m4)

plot(model4)

compareML(model3, model4)

## model3: weight_kg ~ is_boy + age_years + s(ID, bs = "re") + s(ID, age_years, 
## 
## model4: weight_kg ~ is_boy + age_years + s(ID, bs = "re") + s(ID, age_years, 
##  model3:     bs = "re")
## 
## model4:     bs = "re") + s(age_years, ID, bs = "fs", m = 1)
## 
## Chi-square test of fREML scores
## -----
##        Model    Score Edf    Chisq    Df  p.value Sig.
##       model3 9175.449   5                             
## fREML model4 7681.700   7 1493.749 2.000  < 2e-16  ***
## 
## AIC difference: 4720.48, model model4 has lower AIC.

Checking Auto Correlation Function

par(mfrow=c(1,3), cex=1.1)
library(itsadug)
acf_resid(model2, split_pred="ID", main="ACF resid(model2)")
acf_resid(model3, split_pred="ID", main="ACF resid(model3)")
acf_resid(model4, split_pred="ID", main="ACF resid(model4)")

#Weight Prediction using the model 
a_unique_07$pred_wt<-predict(model4,a_unique_07, type='response')

#Comparing Actual and Predicted weight
ggplot(a_unique_07,aes(x=measurement_date, y=pred_wt, group=ID))+
  geom_line()+
  geom_point(aes(x=measurement_date, y=weight_kg, group=ID))

#Checking Actual and predicted for boy and girl
#BOY
boy_wt<-a_unique_07 %>% filter((gender) %in% c("boy"))

weight_mean_boy<-boy_wt %>% group_by(measurement_date) %>% 
  dplyr::summarise(wt_mean=mean(weight_kg),pred_wt_mean=mean(pred_wt))

boy<-ggplot(weight_mean_boy)+
  geom_line(aes(x=measurement_date,y=wt_mean,color='blue'))+
  geom_line(aes(x=measurement_date,y=pred_wt_mean,color='red'))+
  scale_color_discrete(name = "Legends", labels = c( "Predicted","Actual"))

#GIRL

girl_wt<-a_unique_07 %>% filter((gender) %in% c("girl"))

weight_mean_girl<-girl_wt %>% group_by(measurement_date) %>% 
  dplyr::summarise(wt_mean=mean(weight_kg),pred_wt_mean=mean(pred_wt))

girl<-ggplot(weight_mean_girl)+
  geom_line(aes(x=measurement_date,y=wt_mean,color='blue'))+
  geom_line(aes(x=measurement_date,y=pred_wt_mean,color='red'))+
  scale_color_discrete(name = "Legends", labels = c( "Predicted","Actual"))

#Checking mean of actual and predicted weight
weight_mean<-a_unique_07 %>% group_by(measurement_date) %>% 
  dplyr::summarise(wt_mean=mean(weight_kg),pred_wt_mean=mean(pred_wt))

all<-ggplot(weight_mean)+
  geom_line(aes(x=measurement_date,y=wt_mean,color='blue'))+
  geom_line(aes(x=measurement_date,y=pred_wt_mean,color='red'))+
  scale_color_discrete(name = "Legends", labels = c( "Predicted","Actual"))

#Comparing
par(mfrow=c(1,3), cex=.1)
plot_grid(boy, girl, all, labels = c("Weight of Boys","Weight of girl", "Overall weight curve"))

BMI Analysis

#Boys BMI
ggplot(a_unique_07_boys, aes(x=age_years, y=BMI, color=ID)) + geom_line()

#Girls BMI
ggplot(a_unique_07_girls, aes(x=age_years, y=BMI, color=ID)) + geom_line()

model1_bmi <- bam(BMI ~ is_boy + age_years,
          data=a_unique_07)
summary(model1_bmi)

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## BMI ~ is_boy + age_years
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 12.30117    0.29306  41.975  < 2e-16 ***
## is_boy1     -0.57784    0.14132  -4.089 4.45e-05 ***
## age_years    0.63989    0.02184  29.297  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## R-sq.(adj) =  0.224   Deviance explained = 22.5%
## -REML = 8317.3  Scale est. = 14.888    n = 3002

gam.check(model1_bmi)

## 
## Method: REML   Optimizer: outer newton
## full convergence after 5 iterations.
## Gradient range [-2.975128e-07,-2.975128e-07]
## (score 8317.347 & scale 14.88783).
## Hessian positive definite, eigenvalue range [1499.5,1499.5].
## Model rank =  3 / 3

Random intercepts Effect

model2_bmi <- bam(BMI ~ is_boy + age_years
          + s(ID, bs="re"),
          data=a_unique_07)
#summary(m2)

plot(model2_bmi)

gam.check(model2_bmi)

## 
## Method: fREML   Optimizer: perf newton
## full convergence after 8 iterations.
## Gradient range [-2.151038e-06,2.085898e-06]
## (score 6055.329 & scale 2.659805).
## Hessian positive definite, eigenvalue range [64.8164,1502.8].
## Model rank =  144 / 144 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##        k' edf k-index p-value
## s(ID) 141 138      NA      NA

Random intercepts + slopes Effect

model3_bmi <- bam(BMI ~ is_boy + age_years
          + s(ID, bs="re")
          + s(ID, age_years, bs="re"),
          data=a_unique_07)
#summary(m2)

plot(model3_bmi)

gam.check(model3_bmi)

## 
## Method: fREML   Optimizer: perf newton
## full convergence after 9 iterations.
## Gradient range [-4.580215e-09,3.727607e-09]
## (score 5398.14 & scale 1.431173).
## Hessian positive definite, eigenvalue range [51.28346,1505.492].
## Model rank =  285 / 285 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##                  k' edf k-index p-value
## s(ID)           141 131      NA      NA
## s(ID,age_years) 141 132      NA      NA

Comparing Models

library(itsadug)
compareML(model2_bmi,model3_bmi)

## model2_bmi: BMI ~ is_boy + age_years + s(ID, bs = "re")
## 
## model3_bmi: BMI ~ is_boy + age_years + s(ID, bs = "re") + s(ID, age_years, 
##  model2_bmi: BMI ~ is_boy + age_years + s(ID, bs = "re")
## 
## model3_bmi:     bs = "re")
## 
## Chi-square test of fREML scores
## -----
##            Model    Score Edf   Chisq    Df  p.value Sig.
##       model2_bmi 6055.329   4                            
## fREML model3_bmi 5398.140   5 657.189 1.000  < 2e-16  ***
## 
## AIC difference: 1743.39, model model3_bmi has lower AIC.

Random intercepts + slopes + smooths Effect

model4_bmi <- bam(BMI ~ is_boy + age_years
          + s(ID, bs="re")
          + s(ID, age_years, bs="re")
          + s(age_years, ID, bs="fs", m=1),
          data=a_unique_07)
#summary(m4)

plot(model4_bmi)

compareML(model3_bmi, model4_bmi)

## model3_bmi: BMI ~ is_boy + age_years + s(ID, bs = "re") + s(ID, age_years, 
## 
## model4_bmi: BMI ~ is_boy + age_years + s(ID, bs = "re") + s(ID, age_years, 
##  model3_bmi:     bs = "re")
## 
## model4_bmi:     bs = "re") + s(age_years, ID, bs = "fs", m = 1)
## 
## Chi-square test of fREML scores
## -----
##            Model    Score Edf   Chisq    Df  p.value Sig.
##       model3_bmi 5398.140   5                            
## fREML model4_bmi 4671.215   7 726.925 2.000  < 2e-16  ***
## 
## AIC difference: 2515.80, model model4_bmi has lower AIC.

Checking Auto Correlation Function

par(mfrow=c(1,3), cex=1.1)
library(itsadug)
acf_resid(model2_bmi, split_pred="ID", main="ACF resid(model2)")
acf_resid(model3_bmi, split_pred="ID", main="ACF resid(model3)")
acf_resid(model4_bmi, split_pred="ID", main="ACF resid(model4)")

#BMI Prediction using the model 
a_unique_07$pred_bmi<-predict(model4_bmi,a_unique_07, type='response')

#Comparing Actual and Predicted BMI
ggplot(a_unique_07,aes(x=measurement_date, y=pred_bmi, group=ID))+
  geom_line()+
  geom_point(aes(x=measurement_date, y=BMI, group=ID))

#Checking Actual and predicted for boy and girl
#BOY
boy_bmi<-a_unique_07 %>% filter((gender) %in% c("boy"))

bmi_mean_boy<-boy_bmi %>% group_by(measurement_date) %>% 
  dplyr::summarise(bmi_mean=mean(BMI),pred_bmi_mean=mean(pred_bmi))

boy<-ggplot(bmi_mean_boy)+
  geom_line(aes(x=measurement_date,y=bmi_mean,color='blue'))+
  geom_line(aes(x=measurement_date,y=pred_bmi_mean,color='red'))+
  scale_color_discrete(name = "Legends", labels = c( "Predicted","Actual"))

#GIRL

girl_bmi<-a_unique_07 %>% filter((gender) %in% c("girl"))

bmi_mean_girl<-girl_bmi %>% group_by(measurement_date) %>% 
  dplyr::summarise(bmi_mean=mean(BMI),pred_bmi_mean=mean(pred_bmi))

girl<-ggplot(bmi_mean_girl)+
  geom_line(aes(x=measurement_date,y=bmi_mean,color='blue'))+
  geom_line(aes(x=measurement_date,y=pred_bmi_mean,color='red'))+
  scale_color_discrete(name = "Legends", labels = c( "Predicted","Actual"))

#Checking mean of actual and predicted bmi
bmi_mean<-a_unique_07 %>% group_by(measurement_date) %>% 
  dplyr::summarise(bmi_mean=mean(BMI),pred_bmi_mean=mean(pred_bmi))

all<-ggplot(bmi_mean)+
  geom_line(aes(x=measurement_date,y=bmi_mean,color='blue'))+
  geom_line(aes(x=measurement_date,y=pred_bmi_mean,color='red'))+
  scale_color_discrete(name = "Legends", labels = c( "Predicted","Actual"))

par(mfrow=c(1,3), cex=.1)
plot_grid(boy, girl, all, labels = c("BMI of Boys","BMI of girl", "Overall BMI curve"))

#And Compare both shortest and tallest
#check weight and BMI and heartrate as well


#Clustering of growths.
#Model Explaination if not used in the lectures.
#Weight analysis